function filename = calpsd4manual2auto4main(path, file, ouRoot)
    sfrqs = 0.1;
    efrqs = 0.5;
    sfullfile = fullfile(path,file);
    ST = load(sfullfile);
    [pathstr, name, suffix]=fileparts(sfullfile);
    filename = sprintf('%s-SFFT.mat', name);
    %pathname = char('.\TMP\');
    pathname = ouRoot;
    ouFile1 = sprintf('%s%s-Ef.txt',pathname,filename);
    ouFile2 = sprintf('%s%s-Ts.txt',pathname,filename);
    if ~exist(ouFile2) && ~exist(ouFile2)
        if isfield(ST,'FULLDT')
           copypreprocessdata(sfullfile);
           [DATE, VALU, iy, sdate, edate] = readpreprocess();
           %[DATE, VALU, iy, sdate, edate] = readpreprocess4main(sdate, edate);
           winds    = 1024;
           [freq, Ef, spectre0] = connectedDaysPSD( DATE, VALU, winds);
           [TDATE, TDATA] = caltimeserialPSD4NaN(freq, sfrqs, efrqs, spectre0, sdate, edate, winds);
           ouFile = char('RESULT.mat');
           strcmd = char('del RESULT.mat');
           dos(strcmd);
           save(ouFile, 'freq', 'Ef', 'spectre0', 'TDATE', 'TDATA', 'sfrqs', 'efrqs');
           %然后保存
           [pathstr, name, suffix]=fileparts(sfullfile);
           filename = sprintf('%s-SFFT.mat', name);
           %pathname = char('..\TMP\');
           pathname = ouRoot;
           if ~exist(pathname)
              mkdir(pathname);
           end
           strcmd = sprintf('copy RESULT.mat %s%s', pathname, filename);
           dos(strcmd);

           ouFile1 = sprintf('%s%s-Ef.txt',pathname,filename);
           ouFile2 = sprintf('%s%s-Ts.txt',pathname,filename);
           ouputEf2Ts('RESULT.mat', ouFile1,ouFile2);
           %然后保存
        end
    end
end
function copypreprocessdata(sfullfile)
%UNTITLED 此处显示有关此函数的摘要
%   此处显示详细说明
strcmd = sprintf('del %s','TEMP_MONTH.mat');
dos(strcmd);
strcmd = sprintf('copy %s %s',sfullfile,'TEMP_MONTH.mat');
dos(strcmd);
end
%因为观测数据中可能出现缺记导致数据不完整; %这个算法并不好，因为发现很慢需要重新换算法;
function DATA = recovermatdata2index(FULLDT, sdate, edate)
   M = (edate - sdate +1) * 86400;
   [m, n] = size(FULLDT);
   DATA = FULLDT;
   if M ~= m
      DATA = NaN(M, 7);
      AL_M = [1 : M]';
      DATE   = datenum(FULLDT(:,1),FULLDT(:,2),FULLDT(:,3),FULLDT(:,4),FULLDT(:,5),FULLDT(:,6));
      AL_N   = (DATE - DATE(1)) * 86400 + 1; 
      AL_N   = round(AL_N);
      DATA(AL_N,1:7) = FULLDT(:, 1:7);
      ind = find(isnan(DATA(:,1)));
      cdate = sdate + (ind -1)/86400;
      for i = 1 : length(cdate)
          [iy,im,id,iH,iM,iS] = datevec(cdate(i));
          index = ind(i);
          DATA(index,1) = iy;
          DATA(index,2) = im;
          DATA(index,3) = id;
          DATA(index,4) = iH;
          DATA(index,5) = iM;
          DATA(index,6) = iS;
      end
      
%       D_AL_N = diff(AL_N);
%       ind    = find(D_AL_N > 1);
%       ind    = ind + 1;
      
   end
end
function DATA = recovermatdata(FULLDT, sdate, edate)
   M = (edate - sdate +1) * 86400;
   [m, n] = size(FULLDT);
   DATA = FULLDT;
   if M ~= m
       DATA = NaN(M, 8);
       index = 0;
       for i = sdate : edate
           [iy,im,id] = datevec(i);
           for j = 0 : 23
               for k = 0 : 59
                   for l = 0 : 59
                       ind = find(FULLDT(:, 1)==iy & FULLDT(:, 2)==im  & FULLDT(:, 3)==id & ...
                                  FULLDT(:, 4)==j  & FULLDT(:, 5)==k   & FULLDT(:, 6)==l);
                       index = index + 1;
                       DATA(index, 1) = iy;
                       DATA(index, 2) = im;
                       DATA(index, 3) = id;
                       DATA(index, 4) =  j;
                       DATA(index, 5) =  k;
                       DATA(index, 6) =  l;
                       if length(ind) == 1
                          DATA(index, 7) = FULLDT(ind(1), 7);
                          DATA(index, 8) = FULLDT(ind(1), 8);
                       end
                   end
               end
           end
       end       
   end
end

function [DATE, VALU, iy, sdate, edate] = readpreprocess4main(m_sdate, m_edate)
%UNTITLED5 此处显示有关此函数的摘要
%   此处显示详细说明
ST = load('TEMP_MONTH.mat');
DATA = ST.FULLDT;
%DATA = recovermatdata(DATA, m_sdate, m_edate);%这个算法并不好，因为发现很慢需要重新换算法;
DATA = recovermatdata2index(DATA, m_sdate, m_edate);
DATE = datenum(DATA(:,1),DATA(:,2),DATA(:,3),DATA(:,4),DATA(:,5),DATA(:,6));
ind  = find(DATA == 999999);
DATA(ind) = NaN;
index  = firstDATA(DATA(:, 7));
F_VALU = DATA(index, 7);
VALU = DATA(:,7) - F_VALU;
     
sdate = DATE(1);
edate = DATE(end);
[iy,im,id] = datevec(sdate);
[jy,jm,jd] = datevec(edate);
sdate = datenum(iy,im, id);
edate = datenum(jy,jm, jd);
end
function [DATE, VALU, iy, sdate, edate] = readpreprocess()
%UNTITLED5 此处显示有关此函数的摘要
%   此处显示详细说明
ST = load('TEMP_MONTH.mat');
DATA = ST.FULLDT;
DATE = datenum(DATA(:,1),DATA(:,2),DATA(:,3),DATA(:,4),DATA(:,5),DATA(:,6));
ind  = find(DATA == 999999);
DATA(ind) = NaN;
index  = firstDATA(DATA(:, 7));
F_VALU = DATA(index, 7);
VALU = DATA(:,7) - F_VALU;
     
sdate = DATE(1);
edate = DATE(end);
[iy,im,id] = datevec(sdate);
[jy,jm,jd] = datevec(edate);
sdate = datenum(iy,im, id);
edate = datenum(jy,jm, jd);
end
function [freq, Ef, sp0] = connectedDaysPSD( DATE, VALU, winds)
%UNTITLED 此处显示有关此函数的摘要
%   此处显示详细说明
freq     = [];
Ef       = [];
sp0      = [];
N        = length(DATE);
DAYS     = N / 86400;
for i = 1 : DAYS
    s_ind = 1 + (i - 1)* 86400;
    e_ind = i * 86400;
    S_DATE= DATE(s_ind:e_ind);
    S_VALU= VALU(s_ind:e_ind);
   [c_freq, c_Ef, c_sp0] = calmeanPSD4NaN(S_DATE, S_VALU, winds);
    freq = c_freq;
    sp0  = [sp0  c_sp0];
end

ind = find(~isnan(sp0(1, :)));
Ef  = mean(sp0(:, ind), 2);           %计算单边谱
end
function [freq, Ef, sp0] = calmeanPSD4NaN(SEL_SUB_DATE, SEL_SUB_VALU, winds)
 elevation = SEL_SUB_VALU * 10^-8; % 微伽 -> m^2/Hz
 elevation1= elevation;
 nmes      = length(elevation);
 % Below is 3-point filter with weights [ 0.25 0.5 0.25 ] % 4点平均
 elevation(2:nmes-1) = 0.25.*((elevation1(3:nmes)+elevation1(1:nmes-2)) + 2.*elevation1(2:nmes-1));
 Fs  = 1.;                       % sampling frequency in Herz ... should be equal to 1/(times(2)-times(1))
 nfft= winds;                     % nfft specifies the FFT length 
 df  = Fs/nfft;                  % Frequential resolution 
 numoverlap = nfft/2 ;          % numoverlap is the number of samples by which the sections overlap.
 Nf  = nfft/2+1;                   
 Nfo = (Nf-1)/2;
   
 hanning=transpose(0.5 * (1-cos(2 * pi * linspace(0,nfft-1,nfft)/(nfft-1))));  %Vector of Hanning window
 wc2  =1/mean(hanning.^2);                                                   % window correction factor
 freq1=linspace(0,df*(nfft-1),nfft);
 freq =linspace(df,df*(Nf-1),Nf-1)';  
 
 NS1= floor(nmes/nfft);
 NS = NS1*2 - 1;
 Eh = ones(1, NS);
 H  = hanning*Eh;               % array with Hanning windows on each column
 elevmat = zeros(nfft,NS);     % creates array 
 elevmat(1:nfft,1:NS1) = reshape(elevation(1:nfft*NS1),[nfft,NS1]);
 elevmat(1:nfft,NS1+1:NS) = reshape(elevation(1+nfft/2:nfft/2+nfft*(NS1-1)),[nfft,NS1-1]);
 elevmat2=detrend(elevmat,'linear');          % We remove the linear trend from each column
 %Hs1=4.*sqrt(mean(var(elevmat2)));             % We compute the significant wave height 
 elevmat2=H.*detrend(elevmat2,'constant');    % remove mean and apply window
 %Hs2=4.*sqrt(mean(var(elevmat2)));             % Check the impact on the wave height ... must be corrected
  
  % Computes FFT
 zspec=fft(elevmat2,nfft,1)/nfft;
  % computes Power Spectral Density (PSD) 
 Szz=(abs(zspec).^2)*wc2/df;
 % Gathers the variance on both sides of the Nyquist 
 % -> single-sided spectrum
 % 保存频谱分析结果;
 spectre0    = Szz(2:Nf,:);
 spectre0    = acc2disp4sfft(spectre0, freq);
 sp1         = spectre0( : , 1 : NS1);
 sp2         = spectre0( : , 1 + NS1 : 2 * NS1 - 1);
 sp0         = NaN(nfft/2, NS1);
 sp0(:, 1 : NS1 - 1) =  (sp2(:, 1 : NS1 - 1) + sp1(:, 1 : NS1 - 1)) * 0.5;
 sp0(:, NS1)         =   sp1(:, NS1);
 ind = find(~isnan(spectre0(1, :)));
 Ef  = mean(spectre0(:, ind), 2);           %计算单边谱
 Ef  = acc2disp4sfft(Ef, freq);
end
function ouputEf2Ts(inFile, ouFile1,ouFile2)
    %UNTITLED 此处显示有关此函数的摘要
    %   此处显示详细说明
    TS = load(inFile);
    Ef  = TS.Ef;
    freq= TS.freq;
    tdate = TS.TDATE;
    sdate = datestr(tdate, 'yyyy-mm-dd HH:MM:SS');
    tdata = TS.TDATA;

    fid = fopen(ouFile1, 'w');
    line = sprintf('frequency(Hz), Psd(um^2/Hz)');
    fprintf(fid,'%s\r\n',line);
    for i = 1 : length(freq)
       line = sprintf('%10.6f, %10.6f',freq(i), Ef(i) * 2 * 10^12);
       fprintf(fid,'%s\r\n',line);
    end
    fclose(fid);%输出是双边谱
    
    fid = fopen(ouFile2, 'w');
    line = sprintf('YYYY-MM-DD HH:mm:ss, Disp(um)');
    fprintf(fid,'%s\r\n',line);
    for i = 1 : length(sdate)
       line = sprintf('%s, %10.6f', sdate(i,:),tdata(i) * sqrt(2)*10^6);
       fprintf(fid,'%s\r\n',line);
    end
    fclose(fid);

end
function disp4f = acc2disp4sfft(acc4f,frq0)
    i_spectre_log = 10*log10(acc4f);
    [mm,nn]       = size(i_spectre_log);
    period        = 1./frq0;
    transacc2disp = repmat(10 * log10((period/(2 * pi)).^2),[1,nn]);
    i_spectre_p   = i_spectre_log + transacc2disp;
    disp4f        = power(10,i_spectre_p ./ 10.0); % 变为米^2/Hz
end
function [DATE, DATA] = caltimeserialPSD4NaN(freq, sfrqs, efrqs, spectre0, sdate, edate, winds)
%UNTITLED15 此处显示有关此函数的摘要
%   此处显示详细说明
ind    = find(freq >= sfrqs & freq <= efrqs);
dfs    = freq(2) - freq(1);
[m, n] = size(spectre0(ind, :)); 
dfs4a  = repmat(dfs, m, n);

P      = spectre0(ind, :) .* dfs4a;
DATA   = sqrt(sum(P)) ;
DATE   = calSplit4Days(sdate, edate, winds);
% 计算的平均功率谱 双边谱;
end

function index = firstDATA(VALU)
   index = -1;
   N = length(VALU);
   for i = 1 : N
       if ~isnan(VALU(i))
          index = i;
          break;
       end
   end
end

function DATE4SERIAL = calSplit4Days(sdate, edate, winds)
%    SDATE = datenum(iy,im,1);
%    EDATE = addtodate(SDATE,1,'month');
%    EDATE = addtodate(EDATE,-1,'day');
   N     = edate - sdate + 1;
   N4day = floor(86400 / winds);
   NUM   = N4day * N;
   DATE4SERIAL = NaN(NUM, 1);
   index       = 1;
   for i = 1 : N
       TDATE = addtodate(sdate, i - 1, 'day');
       TDATE = addtodate(TDATE, winds/2, 'second');
       for j = 1 : N4day
           CDATE = addtodate(TDATE, (j - 1) * winds, 'second');
           DATE4SERIAL(index) = CDATE;
           index = index + 1;
           if index > NUM
              break;
           end
       end
   end
end